---
title: "Political Responsibility in the Era of COVID-19" 
author: "Haruka Braun, Lucy Ding, and Camryn Jones"
date: '\today'
fontsize: 12pt
output:
  pdf_document: default
urlcolor: red
header-includes:
  - \usepackage{bm}
  - \usepackage{placeins}
  - \usepackage[labelformat=empty]{caption}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r, warning = FALSE, message = FALSE}
# Load necessary packages
library(usmap)
library(haven)
library(tidyverse)
library(dplyr)
library(sf)
library(ggplot2)
library(janitor)
library(sandwich)
library(lmtest)
library(AER)
library(stargazer)
library(lfe)
library(estimatr)
library(stats)
library(MatchIt)
library(optmatch)
library(tableone)
library(cobalt)
library(CBPS)
library(gt)
library(tigris)
library(gtsummary)
library(psych)
library(fixest)
library(sjPlot)
library(sjmisc)
library(sjlabelled)
library(modelsummary)
options(scipen = 100)
select <- dplyr::select
knitr::opts_chunk$set(tidy.opts=list(width.cutoff=80), tidy=TRUE)
```

```{r, warning = FALSE, message = FALSE}
# Read in paper data and election data
Data <- read_dta("repli_main_tables.dta") %>%
  clean_names()
Election_2020 <- read_csv("2020_0_0_2.csv") %>%
  clean_names()
Election_2016 <- read_csv("county_2016.csv") %>%
  clean_names()

# Merge paper data with election data
Final <- merge(Data, Election_2020, by = "fips")
Final <- merge(Final, Election_2016, by = "fips")

# Load fips data and merge with final data
data(fips_codes)
fips_codes <- fips_codes %>%
  mutate(fips = paste(state_code, county_code, sep = "")) %>%
  mutate(fips = as.numeric(fips))
Final <- left_join(Final, fips_codes, by = "fips")

# Clean final data
Final <- Final %>%
  clean_names() %>%
  mutate(
    total_vote_20 = as.numeric(total_vote_x),
    trump_20 = as.numeric(donald_j_trump_x),
    total_vote_16 = as.numeric(total_vote_y),
    trump_16 = as.numeric(donald_j_trump_y),
    clinton_16 = as.numeric(hillary_clinton),
    biden = as.numeric(joseph_r_biden_jr)
  ) %>%
  mutate(state_code_num = as.numeric(state_code)) %>%
  filter(state_code_num < 60) %>%
  mutate(
    trump_voting_2020 = (trump_20 / total_vote_20) * 100,
    changes_trump = (trump_20 / total_vote_20 - trump_16 / total_vote_16) * 100,
    changes_total = total_vote_20 - total_vote_16,
    covid_cases = (population / 10000) * cases_10000,
    covid_deaths = (population / 100000) * deaths_100000,
    rep_share_2020 = (trump_20 / total_vote_20) * 100,
    rep_share_2016 = (trump_16 / total_vote_16) * 100,
    depvar = rep_share_2020 - rep_share_2016,
    covid_cases = (population / 10000) * cases_10000,
    cases_10000_2 = cases_10000 * cases_10000,
    trump_gap = trump_20 - biden
  )
```

```{r, message = FALSE}
#Figure 1
plot_usmap(data = Final, values = 'cases_oct_10000', exclude = c('AK', 'HI'),
           size = 0.1) + 
  scale_fill_continuous() + scale_fill_gradientn(colors = c("white", 
                                                            "dodgerblue4", 
                                                            "black"), 
                                                 name = "Cases per 10,000") +
  theme(legend.position = "left") + 
  labs(title = "Fig 1. Cumulative Number of COVID-19 cases per 10,000")
```

```{r, warning = FALSE, message = FALSE}
#Figure 2
plot_usmap(
  data = Final, values = "depvar", exclude = c("AK", "HI"),
  size = 0.1) + scale_fill_continuous() +
  scale_fill_gradientn(colors = c("white", "dodgerblue4", "black"), name =
                         "Change in Vote Share") + 
  theme(legend.position = "left") +
  labs(title =
      "Fig 2. Changes in share of votes for Donald Trump from 2016 to 2020"
  )
```

```{r, warning = FALSE, message = FALSE}
# Drop missing values
tab1df <- Final %>%
  drop_na(trump_voting_2020, workplaces_april, s_emp_s2018) %>%
  select(
    trump_voting_2020, changes_trump, changes_total, covid_cases,
    cases_10000, covid_deaths, deaths_100000, i_animal_processing,
    unemp_change_covid1
  )
```

```{r, warning = FALSE, message = FALSE}
# Create summary table of election outcomes, COVID-19, and labor outcomes
describe(tab1df, fast = TRUE, skew = FALSE) %>%
  select(mean, sd, max, min, n) %>%
  gt(rownames_to_stub = TRUE) %>%
  tab_header(title = "Table 1") %>%
  tab_row_group(
    label = "Labor outcomes",
    rows = c("i_animal_processing", "unemp_change_covid1")
  ) %>%
  tab_row_group(
    label = "COVID-19 incidence",
    rows = c("covid_cases", "cases_10000", "covid_deaths", "deaths_100000")
  ) %>%
  tab_row_group(
    label = "Election outcomes",
    rows = c("trump_voting_2020", "changes_trump", "changes_total")
  ) %>%
  cols_label(
    mean = "Mean",
    sd = "S.D.",
    max = "Max",
    min = "Min"
  ) %>%
  fmt_number(
    columns = 1:5,
    decimals = 5
  )
```

```{r, warning = FALSE, message = FALSE}
# Drop missing values
Final_tab2 <- Final %>%
  drop_na(workplaces_april, s_emp_s2018)

# Column 1 replication - OLS fixed-effects with demographic and socioeconomic
# controls
model_c1 <- felm(data = Final_tab2, depvar ~ cases_10000 + black2018 +
  white_nonhispanic2018 + college2018 + foreignborn2018 +
  female2018 + totalpop2018 + age_20_242010 + age_25_342010 +
  age_35_442010 + age_45_542010 + age_55_592010 +
  age_60_642010 + age_65_742010 + age_75_842010 +
  age_85over2010 + exposure + proximity + remindex +
  crit_work | istate | 0 | istate)

# Column 2 replication - OLS fixed-effects with demographic, socioeconomic,
# and social mobility controls
model_c2 <- felm(data = Final_tab2, depvar ~ cases_10000 + black2018 +
  white_nonhispanic2018 + college2018 + foreignborn2018 +
  female2018 + totalpop2018 + age_20_242010 + age_25_342010 +
  age_35_442010 + age_45_542010 + age_55_592010 +
  age_60_642010 + age_65_742010 + age_75_842010 +
  age_85over2010 + exposure + proximity + remindex +
  crit_work + workplaces_april | istate | 0 | istate)

# Column 3 replication - OLS fixed-effects with demographic, socioeconomic,
# social mobility, and unemployment controls
model_c3 <- felm(data = Final_tab2, depvar ~ cases_10000 + black2018 +
  white_nonhispanic2018 + college2018 + foreignborn2018 +
  female2018 + totalpop2018 + age_20_242010 + age_25_342010 +
  age_35_442010 + age_45_542010 + age_55_592010 +
  age_60_642010 + age_65_742010 + age_75_842010 +
  age_85over2010 + exposure + proximity + remindex +
  crit_work + workplaces_april + unemp_change_covid1 |
  istate | 0 | istate)
```

```{r, warning = FALSE, message = FALSE}
# Table 2 replication with standard errors and 90% confidence intervals
models <- list(model_c1, model_c2, model_c3)
display <- c("cases_10000"  = "Cumulative COVID Cases", "exposure" = "Exposure", "proximity" = "Proximity", "remindex" = "Remote Work", "crit_work" = "Essential Work", "workplaces_april" = "April 1st Acvitivity Change", "unemp_change_covid1" = "Unemployment")
demorow <- tribble(~term, ~model_c1, ~model_c2, ~model_c3,
                   'Demographic Controls', 'Yes', 'Yes', 'Yes',
                   'Socioeconomic Controls', 'Yes', 'Yes', 'Yes')
attr(demorow, 'position') <- c(15, 16)
table2 <- modelsummary(models,
             title = "Table 2 Replication with Standard Errors and 90% Confidence Intervals",
             fmt= 5,
             output = 'gt',
             stars = TRUE,
             estimate = "{estimate} ({std.error})",
             statistic = "[{conf.low}, {conf.high}]",
             coef_map = display,
             add_rows = demorow) %>% tab_header(title = "Table 2")
table2
```

```{r, warning = FALSE, message = FALSE}
# Table 3 Replication: Counterfactual Analysis with 90% confidence intervals
c1 <- confint(model_c1, level = 0.90)
Final %>%
  mutate(paper = -0.0012 * cases_10000 * -0.05,
         c_upper = c1[1] * cases_10000 * -0.05,
         c_lower = c1[1,2] * cases_10000 * -0.05,
         additional = paper * total_vote_20,
         additional_upper = c_upper * total_vote_20,
         additional_lower = c_lower * total_vote_20) %>%
  group_by(state) %>%
  summarize(trump_gap = sum(trump_gap),
            additional = sum(additional),
            upper_bound = sum(additional_upper),
            lower_bound = sum(additional_lower)
            ) %>%
  filter(state %in% c("AZ", "GA", "MI", "PA", "WI")) %>%
  gt() %>% 
  tab_header(title = "Table 3") %>%
  cols_label(
    state = md("State"),
    trump_gap = md("Trump Gap"),
    additional = md("5%"),
    upper_bound = md("CI Upper Bound"),
    lower_bound = md("CI Lower Bound")
  )
```

```{r, warning = FALSE, message = FALSE}
# Read in county level data for the US House of Representatives
county_house_2020 <- read_csv("county_house_2020.csv") %>%
  clean_names()
county_house_2018 <- read_csv("county_house_2018.csv") %>%
  clean_names()

# Join with control variables and clean data
county_house <-
  left_join(county_house_2020, county_house_2018,
    by = c("fips", "geographic_name", "geographic_subtype")
  ) %>%
  slice(-c(1)) %>%
  mutate(
    fips = as.numeric(fips),
    total_vote_20 = as.numeric(total_vote.x),
    democratic_20 = as.numeric(democratic.x),
    republican_20 = as.numeric(republican.x),
    total_vote_18 = as.numeric(total_vote.y),
    democratic_18 = as.numeric(democratic.y),
    republican_18 = as.numeric(republican.y)
  ) %>%
  select(
    fips, geographic_name, total_vote_20, democratic_20, republican_20,
    total_vote_18, democratic_18, republican_18
  )

extension <- left_join(Data, county_house, by = "fips") %>%
  mutate(
    rep_vote_share_20 = (republican_20 / total_vote_20) * 100,
    rep_vote_share_18 = (republican_18 / total_vote_18) * 100,
    depvar_ext = rep_vote_share_20 - rep_vote_share_18
  )
```
\newpage
```{r, message = FALSE}
# Fixed effects linear model identical to column 1 of Table 2 using data from the US House
model_ext <- felm(data = extension %>% drop_na(workplaces_april, s_emp_s2018), 
                  depvar_ext ~ cases_10000 + black2018 + white_nonhispanic2018 + 
                    college2018 + foreignborn2018 + female2018 + totalpop2018 + 
                    age_20_242010 + age_25_342010 + age_35_442010 + 
                    age_45_542010 + age_55_592010 + age_60_642010 + 
                    age_65_742010 + age_75_842010 + age_85over2010 + exposure + 
                    proximity + remindex + crit_work | istate | 0 | istate)

model_extm <- list("Model" = model_ext)
display2 <- c("cases_10000"  = "Cumulative COVID Cases", "exposure" = "Exposure", "proximity" = "Proximity", "remindex" = "Remote Work", "crit_work" = "Essential Work")
demrow <- tribble(~term, ~model_c1, 
                   'Demographic Controls', 'Yes',
                  'Socioeconomic Controls', 'Yes')
attr(demrow, 'position') <- c(15,16)
table4 <- modelsummary(model_extm,
             output = 'gt',
             fmt= 5,
             coef_map = display2,
             add_rows = demrow) %>% tab_header(title = "Table 4")
table4

```


```{r, warning = FALSE, message = FALSE}
# Read in county level data for Governors
county_governors_2016 <- read_csv("2016 governor.csv") %>%
  clean_names()
county_governors_2020 <- read_csv("2020 governor.csv") %>%
  clean_names()

# Join with control variables and clean data
county_governors <-
  left_join(county_governors_2020, county_governors_2016,
    by = c("fips", "geographic_name", "geographic_subtype")
  ) %>%
  slice(-c(1)) %>%
  mutate(
    fips = as.numeric(fips),
    total_vote_20 = as.numeric(total_vote.x),
    democratic_20 = as.numeric(democratic.x),
    republican_20 = as.numeric(republican.x),
    total_vote_16 = as.numeric(total_vote.y),
    democratic_16 = as.numeric(democratic.y),
    republican_16 = as.numeric(republican.y)
  ) %>%
  select(
    fips, geographic_name, total_vote_20, democratic_20, republican_20,
    total_vote_16, democratic_16, republican_16
  )

gextension <- left_join(Data, county_governors, by = "fips") %>%
  mutate(
    rep_vote_share_20 = (republican_20 / total_vote_20) * 100,
    rep_vote_share_16 = (republican_16 / total_vote_16) * 100,
    depvar_ext = rep_vote_share_20 - rep_vote_share_16,
    rep_gap = republican_20 - democratic_20
  )

# Fixed effects linear model identical to column 1 of Table 2 using data from the Election of Governors
model_ext <- felm(data = gextension %>% drop_na(workplaces_april, s_emp_s2018), 
                  depvar_ext ~ cases_10000 + black2018 + white_nonhispanic2018 + 
                    college2018 + foreignborn2018 + female2018 + totalpop2018 + 
                    age_20_242010 + age_25_342010 + age_35_442010 + 
                    age_45_542010 + age_55_592010 + age_60_642010 + 
                    age_65_742010 + age_75_842010 + age_85over2010 + exposure + 
                    proximity + remindex + crit_work  | istate | 0 | istate)


model_extm <- list("Model" = model_ext)
demrow <- tribble(~term, ~model_c3, 
                   'Demographic Controls', 'Yes',
                  'Socioeconomic Controls', 'Yes')
attr(demrow, 'position') <- c(15,16)
table5<- modelsummary(model_extm,
             fmt= 5,
             output = 'gt',
             coef_map = display2,
             add_rows = demrow) %>% tab_header(title = "Table 5")
table5
```

```{r, warning = FALSE, message = FALSE}
# Counterfactuals with 90% confidence intervals
statecountyfips <- read.csv('fips_by_state.csv')
gextension <- merge(gextension, statecountyfips, by='fips')
#c11 <- confint(model_ext, level = 0.90)
gextension %>%
  mutate(paper1 = -0.00207 * cases_10000 * -0.05,
         paper2 = -0.00207 * cases_10000 * -0.1,
         paper3 = -0.00207 * cases_10000 * -0.2,
         additional_05 = paper1 * total_vote_20,
         additional_10 = paper2 * total_vote_20,
         additional_20 = paper3 * total_vote_20) %>%
  group_by(state) %>%
  summarize(rep_gap = sum(rep_gap),
            additional_05 = sum(additional_05),
            additional_10 = sum(additional_10),
            additional_20 = sum(additional_20)
            )%>%
  filter(state %in% c("DE", "IN", "MO", "MT", "NC", "ND", "NH", "UT",
                      "VT", "WA", "WV")) %>% 
  gt() %>% tab_header(title = "Table 6")%>%
  cols_label(
    state = md("State"),
    rep_gap = md("Republican Gap"),
    additional_05 = md("5%"),
    additional_10 = md("10%"),
    additional_20 = md("20%")
  )
```